function HJs = Js(Nx,Ny,J1,J2)
%JS Generate the hopping term of Hamiltonian with periodic boundary
%condition
N1=Nx*Ny; % on-lattice
N2=(Nx-1)*(Ny-1); % off-lattice
N=N1+N2;
HJs=zeros(N);

for ii=1:Ny
    if ii~=Ny
        if mod(ii,2)==1 % odd
            for jj=1:Nx-1
                [HJs((ii-1)*Nx+jj,(ii-1)*Nx+jj+1),...
                 HJs((ii-1)*Nx+jj+1,(ii-1)*Nx+jj),...
                 HJs((ii-1)*Nx+jj,ii*Nx+jj),...
                 HJs(ii*Nx+jj,(ii-1)*Nx+jj),...
                 HJs((ii-1)*Nx+jj+1,ii*Nx+jj),...
                 HJs(ii*Nx+jj,(ii-1)*Nx+jj+1)]=deal(J1);
                [HJs((ii-1)*Nx+jj,N1+(ii-1)*(Nx-1)+jj),...
                 HJs(N1+(ii-1)*(Nx-1)+jj,(ii-1)*Nx+jj),...
                 HJs((ii-1)*Nx+jj+1,N1+(ii-1)*(Nx-1)+jj),...
                 HJs(N1+(ii-1)*(Nx-1)+jj,(ii-1)*Nx+jj+1),...
                 HJs(ii*Nx+jj,N1+(ii-1)*(Nx-1)+jj),...
                 HJs(N1+(ii-1)*(Nx-1)+jj,ii*Nx+jj)]=deal(J2);
            end
            [HJs(ii*Nx,(ii+1)*Nx),...
             HJs((ii+1)*Nx,ii*Nx),...
             HJs(ii*Nx,(ii-1)*Nx+1),...
             HJs((ii-1)*Nx+1,ii*Nx),...
             HJs((ii+1)*Nx,(ii-1)*Nx+1),...
             HJs((ii-1)*Nx+1,(ii+1)*Nx)]=deal(J1);
        else % even
            for jj=1:Nx-1
                [HJs((ii-1)*Nx+jj,(ii-1)*Nx+jj+1),...
                 HJs((ii-1)*Nx+jj+1,(ii-1)*Nx+jj),...
                 HJs((ii-1)*Nx+jj,ii*Nx+jj+1),...
                 HJs(ii*Nx+jj+1,(ii-1)*Nx+jj),...
                 HJs((ii-1)*Nx+jj+1,ii*Nx+jj+1),...
                 HJs(ii*Nx+jj+1,(ii-1)*Nx+jj+1)]=deal(J1);
                [HJs((ii-1)*Nx+jj,N1+(ii-1)*(Nx-1)+jj),...
                 HJs(N1+(ii-1)*(Nx-1)+jj,(ii-1)*Nx+jj),...
                 HJs((ii-1)*Nx+jj+1,N1+(ii-1)*(Nx-1)+jj),...
                 HJs(N1+(ii-1)*(Nx-1)+jj,(ii-1)*Nx+jj+1),...
                 HJs(ii*Nx+jj+1,N1+(ii-1)*(Nx-1)+jj),...
                 HJs(N1+(ii-1)*(Nx-1)+jj,ii*Nx+jj+1)]=deal(J2);
            end
            [HJs(ii*Nx,(ii+1)*Nx),...
             HJs((ii+1)*Nx,ii*Nx),...
             HJs(ii*Nx,ii*Nx+1),...
             HJs(ii*Nx+1,ii*Nx),...
             HJs((ii+1)*Nx,ii*Nx+1),...
             HJs(ii*Nx+1,(ii+1)*Nx)]=deal(J1);
        end
    elseif ii==Ny
        for jj=1:Nx-1
            [HJs((ii-1)*Nx+jj,(ii-1)*Nx+jj+1),...
             HJs((ii-1)*Nx+jj+1,(ii-1)*Nx+jj),...
             HJs((ii-1)*Nx+jj,jj+1),...
             HJs(jj+1,(ii-1)*Nx+jj),...
             HJs((ii-1)*Nx+jj+1,jj+1),...
             HJs(jj+1,(ii-1)*Nx+jj+1)]=deal(J1);
        end
        [HJs(N1,Nx),...
         HJs(Nx,N1),...
         HJs(N1,1),...
         HJs(1,N1),...
         HJs(Nx,1),...
         HJs(1,Nx)]=deal(J1);        
    end
end


end